function EulerFluxes(Q, γ)
    """
    Purpose: evaluate primitive variables and Euler flux functions
    """
    # extract conserved variables
    ρ  = Q[:,:,1]
    ρu = Q[:,:,2]
    ρv = Q[:,:,3]
    E  = Q[:,:,4]

    # compute primitive variables
    @inbounds u = ρu./ρ
    @inbounds v = ρv./ρ
    @inbounds p = (γ - 1)*(E - 0.5*(ρu.*u + ρv.*v))

    # Compute flux functions
    F = zeros(size(Q))
    @inbounds F[:,:,1] = ρu
    @inbounds F[:,:,2] = ρu.*u + p
    @inbounds F[:,:,3] = ρv.*u
    @inbounds F[:,:,4] = u.*(E + p)

    G = zeros(size(Q))
    @inbounds G[:,:,1] = ρv
    @inbounds G[:,:,2] = ρu.*v
    @inbounds G[:,:,3] = ρv.*v + p
    @inbounds G[:,:,4] = v.*(E + p)

    F, G, ρ, u, v, p
end
